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Abstract — In recent years, two important techniques for geometric numerical discretization 
have been developed. In computational electromagnetics, spatial discretization has been im- 
proved by the use of mixed finite elements and discrete differential forms. Simultaneously, the 
dynamical systems and mechanics communities have developed structure-preserving time inte- 
grators, notably variational integrators that are constructed from a Lagrangian action principle. 
Here, we discuss how to combine these two frameworks to develop variational spacetime integra- 
tors for Maxwell's equations. Extending our previous work, which first introduced this variational 
perspective for Maxwell's equations without sources, we also show here how to incorporate free 
sources of charge and current. 



1. INTRODUCTION 

In computational electromagnetics, as in an increasing number of other fields in applied science 
and engineering, there is both practical and theoretical interest in developing geometric numerical 
integrators. These numerical methods preserve, by construction, various geometric properties and 
invariants of the continuous physical systems that they approximate. This is particularly important 
for applications where even high-order methods may fail to capture important features of the un- 
derlying dynamics. In this short paper, we show that the traditional Yee scheme and extensions can 
be derived from the Euler-Lagrange equations of a discrete action, i.e., by designing an electromag- 
netic variational integrator, including free sources of charge and current in non-dissipative media. 
Furthermore, we present how to use this discrete geometric framework to allow for asynchronous 
time stepping on unstructured grids, as recently introduced in Stern et al. |T0]. 

Variational integrators (not to be confused with variational methods such as finite element 
schemes) were originally developed for geometric time integration, particularly to simulate dynam- 
ical systems in Lagrangian mechanics. The key idea is the following: rather than approximating 
the equations of motion directly, one discretizes the Lagrangian and its associated action integral 
(e.g., using a numerical quadrature rule), and then derives a structure-preserving approximation to 
the equations of motion by applying Hamilton's principle of stationary action. Since the numerical 
method is derived from a Lagrangian variational principle, some important results from Lagrangian 
dynamics carry over to the discretized system, including Noether's theorem relating symmetries 
to conserved momentum maps, as well as the fact that the Euler-Lagrange flow is a symplectic 
mapping. (Marsden and West [SeeE], Lew et al. [See [7].) 

Overview. To develop a variational integrator for Maxwell's equations, the discrete Hamilton's 
principle needs to incorporate more than just the time discretization, as in mechanics; spatial 
discretization also needs to be handled carefully. Building upon mixed finite elements in space |9l 
[21 15] , we treat the electromagnetic Lagrangian density as a discrete differential 4-form in spacetime. 
Extremizing the integral of this Lagrangian density with fixed boundary conditions directly leads 
to discrete update rules for the electromagnetic fields, with either uniform or asynchronous time 
steps across the various spatial elements. 

2. REVIEW OF MAXWELL'S EQUATIONS IN SPACETIME 

Electromagnetic Forms. Let A be a 1-form on spacetime, called the electromagnetic potential, 
and then define the Faraday 2-form to be its exterior derivative F = dA. Given a time coordinate 
t, this splits into the components 

F = E Adt + B, 
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where E is the electric displacement 1-form and B is the magnetic flux 2-form, both defined on 
the spacelike Cauchy surfaces S with constant t. If * is the Hodge star associated to the spacetime 
metric, then we can also split the dual 2-form 

*F = Adt- *eE = H Adt- D, 

where (again, restricted to Cauchy surfaces) H is the magnetic displacement 1-form, D is the electric 
flux 2-form, and *^ and *c arc respectively the magnetic permeability and electric permittivity. 
Finally, for systems with free sources, there is a source 3-form J', satisfying the continuity of 
charge condition dj' = 0. In terms of coordinates, this can be split into 

J = J A dt - p, 

where J is the current density 2-form and p is the charge density 3-form on Cauchy surfaces. 

Maxwell's Equations. With the spacetime forms and operators defined above. Maxwell's equa- 
tions become 

dF = 0, d*F = J. 

Note that the first equation follows automatically from F = dA, since taking the exterior derivative 
of both sides yields dF = ddA = 0. The second equation is consistent with the continuity of charge 
condition, since dJ = dd*F = 0. 

LagrEingian Formulation. Given the electromagnetic potential 1-form A and source 3-form we 
can define the Lagrangian density to be the 4-form 

C = -^dA A *dA + AAJ, 

with the associated action functional S[A] = J-^ C taken over the spacetime domain X. Suppose 
that a is a variation of A, vanishing on the boundary dX. Varying the action along a yields 

dS[A] • a = / {-da A *dA + a A J) = / a A {-d*dA + J) . 
Jx Jx 

Hamilton's principle of stationary action states that this variation must equal zero for any such a, 
implying the Euler-Lagrange equations d*dA = J . Finally, substituting F = dA and recalling that 
dF = ddA = 0, we see that this is equivalent to Maxwell's equations. 

3. GEOMETRIC PROPERTIES OF MAXWELL'S EQUATIONS 

As written in terms of F above, Maxwell's equations have 8 components: 6 dynamical equations, 
which describe how the fields change in time, and 2 "divergence constraints" containing only spatial 
derivatives. The fact that these constraints are automatically preserved by the dynamical equations 
(and can therefore effectively be ignored except at the initial time) comes directly from the differen- 
tial gauge symmetry and Lagrangian variational structure. We discuss these geometric properties 
here, with a view towards developing numerical methods that preserve them. 

Reduction by Gauge Fixing. Maxwell's equations are invariant under gauge transformations 
A A df foi any scalar function /, since taking the exterior derivative maps F h- > F H- ddf = F. 
Therefore, given a time coordinate t, we can fix the gauge so that A ■ ^ =0, i.e., A has only 
spacclikc components. This partial gauge fixing is known as the Weyl gauge. Restricted to this 
subspace of potentials, the Lagrangian then becomes 

£ = {dtA + d-sA) A * {dtA + dj:A) + AAJ 

= {dtA A *dtA + dj:A A *dT.A) + AaJ Adt 

Here, wc have adopted the notation dt and d-^ for the exterior derivative taken only in time and in 
space, respectively; in particular, we then have dtA = E Adt and d^A = B. 
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Next, varying the action along a restricted variation a that vanishes on (9X, 

AS\A\ ■a= {dta AD - dj^a A H A dt + a A J A dt) (1) 
Jx 

= / a A {dtD - d^H A dt + J A dt) . 
Jx 

Setting this equal to zero by Hamilton's principle, one immediately gets Ampere's law as the sole 
Euler-Lagrange equation. The divergence constraint dsD = p, corresponding to Gauss' law, has 
been eliminated via the restriction to the Weyl gauge. 

Noether's Theorem Automatically Preserves Gauss' Law. There are two ways that one can see 
why Gauss' law is automatically preserved, even though it has been eliminated from the Euler- 
Lagrange equations. The first is to take the "divergence" of Ampere's law, obtaining 

= dj^dtD - dj:dY.H Adt + dT.J A dt = dt {dj^D - p) . 

Therefore, if this condition holds at the initial time, then it holds for all time. 

A more "geometric" way to obtain this result is to use Noether's theorem, with respect to the 
remaining gauge symmetry A A + d-^f for scalar functions / on S. To derive this, let us restrict 
A to be an Euler-Lagrange solution in the Weyl gauge, but remove the previous requirement that 
variations a be fixed at the initial time to and final time tj. Then, varying the action along this 
new a, the Euler-Lagrange term disappears, but we now pick up an additional boundary term due 
to integration by parts 

r */ 
dS[A] ■a= I aAD . 

to 

If we vary along a gauge transformation a = dsf, then this becomes 



dS[A] • ds/ = / ds/ A D 



E 



fAd^D 



f 



Alternatively, plugging a = d^f into [Equation 1[ we get 

dS[A] -d^f = [ dj^f AJ Adt = - I f Ad^J Adt = - [ f A dtp = - [ f Ap \ 
Jx Jx Jx Jt, to 

Since these two expressions are equal, and / is an arbitrary function, it follows that 

{d^D-p)\%=0. 

This indicates that dj]D — p is a conserved quantity, a momentum map, so if Gauss' law holds at 
the initial time, then it holds for all subsequent times as well. 

4. GEOMETRIC DISCRETIZATION OF MAXWELL'S EQUATIONS 

Discretizing Maxwell's equations, while preserving the geometric properties mentioned above, can 
be achieved using cochains as discrete substitutes for differential forms, as previously done in, 
e.g., Bossavit [2 . Therefore, to compute Maxwell's equations, we begin by discretizing the 2- 
form F on a spacetime mesh K: F assigns a real value to each oriented 2-face of the mesh. The 
exterior derivative d is discretized by the coboundary operator, so the equation dF = states that 
{dF,a^) = {F,da^) = 0, where is any oriented 3-cell in K and da^ is its 2-chain boundary. 

Next, given a discrete Hodge star operator |3l [TTl [6l [1], *F is a 2-form on the dual mesh *K, 
while is defined as a discrete dual 3-form. Then, for every dual 3-cochain (where is the 
corresponding primal edge), the equation d*F = J becomes 

When the cells and *(T^ are spacetimelike, then these correspond to the dynamical components 
of Maxwell's equations, and can be used to compute subsequent values of F . When the cells are 
purely spacelike, they correspond to the divergence constraint equations. The exact expression and 
update of these fields in time now depends on which type of mesh and time stepping method is 
desired, as described next. 
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Uniform Time Stepping. For uniform rectangular meshes aligned with the {x, y, z, t) axes, we 
can emulate the smooth coordinate expression of F as 

F = Ea:dx Adt + Eydy Adt + E^ dz A dt 

+ Bx dy A dz + By dz Adx + B^ dx A dy. 



This suggests the following discretization of F, shown in [Figure 1 store E^AxAt on the xt-faces 



EyAyAt on the yt-faces, and E^AzAt on the zt-faces; likewise for B. To discretize J', we can 
similarly store J^AyAzAt, JyAzAxAt, J^AxAyAt, and pAxAyAz on the corresponding dual 
3-cells. If we then enforce the equations dF = and d*F = J' in the discrete sense, the result 




Figure 1: The 2-form F = E A dt + B can be discretized, on a rectangular spacetime mesh, by storing the 
components of E and B on 2D faces. The resulting numerical method is Yee's FDTD scheme. 

is precisely the finite-difference time-domain (FDTD) integration scheme of Yee |12] . A similar 
procedure can be applied on unstructured (e.g., simplicial) spatial grids, on which we take uniform 
time steps At (creating prism-shaped spacetime primal elements); in this case, solving the discrete 
Maxwell's equations recovers the more recent "Yee- like" method of Bossavit and Kettunen 



Asynchronous Time Stepping. As initially mentioned in Stern et al. \LQi, a more flexible inte- 
gration scheme can be designed by assigning a different time step per element, in order to focus 
computational power where needed. A visualization of how to store F on such a mesh structure is 



shown in Figure 2] This defines an asynchronous variational integrator that preserves the numerical 
properties of the uniform-stepping methods outlined above. The procedure to repeatedly update 
E and B asynchronously in time is as follows: 

1. Select the face for which B needs to be updated next. 

2. E advances B, using dF = 0. 

3. B advances E on neighboring edges, using d*F = J . 











dx A di 










It 


dx A dt 












E Adt 



Figure 2: Here, F is discretized on an asynchronous time-stepping grid, where each spatial element takes a 
different-sized time step from its neighbors. This can be done for either a rectangular spatial grid (left) or 
an unstructured/simplicial spatial mesh (right). 
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(a) regular grid, uniform time steps 



(b) random grid, asynchronous time steps 



Figure 3: Our geometric integrator robustly maintains near-conservation of energy, even for asynchronous 
time stepping on a random spatial grid. 



Details of this algorithm, along with initial numerical results, can be found in Stern et al. |10] 
where only the case ^7 = was described. 



One promising result from the initial numerical experiments, shown in Figure 3, is that this 
asynchronous integrator matches the FDTD scheme's excellent energy conservation behavior, even 
on a highly irregular grid, without exhibiting artificial damping or forcing. Future work is expected 
to include formal accuracy and stability analysis of this asynchronous integrator, focusing both on 
the choice of Hodge star and on resonance stability criteria for selecting the individual time steps. 
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